Discreteness of populations enervates biodiversity in evolution 
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Abstract 

^-^ Biodiversity widely observed in ecological systems is attributed to the dynamical balance among 

^^ the competing species. The time- varying populations of the interacting species are often captured 
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rather well by a set of deterministic replicator equations in the evolutionary game theory. However, 









O^ 



> 

en 
m 



intrinsic fluctuations arisen from the discreteness of populations lead to stochastic derivations from 
the smooth evolution trajectories. The role of these fluctuations is shown to be critical at causing 
extinction and deteriorating the biodiversity of ecosystem. We use children's rock-paper-scissors 



P^ game to demonstrate how the intrinsic fluctuations arise from the discrete populations and why the 

Q biodiversity of the ecosystem decays exponentially, disregarding the detail parameters for competing 

1^ mechanism and initial distributions. The dissipative trend in biodiversity can be analogized to the 

gradual erosion of kinetic energy of a moving particle due to air drag or fluid viscosity. The 
dissipation-fluctuation theorem in statistical physics seals the fate of these originally conserved 
quantities. This concept in physics can be generalized to scrutinize the errors that might be 
incurred in the ecological, biological, and quantitative economic modeling for which the ingredients 



•^ are all discrete in number. 
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Biodiversity is commonly used to indicate the stability of an ecosystem [TH3]. One of 
the central issues is to effectively promote the biodiversity while attracting more scientists 
attention from various fields[lHH]. The causes that threaten the biodiversity, for instance, 
climate changejUE], over-harvesting, habitat destruction [6], and population mobilityjT], are 
well studied. Above those factors, Darwin's theory of natural selection plays a crucial role 
in catalysis (914T2]. People are warned to reduce these effects in order to maintain and reserve 
the nature's biodiversity. Nevertheless, a naive reversed statement can be checked, namely, 
without any hazardous factors, would ecosystem be perfectly stable? 

Here, we show the emergence of an intrinsic force originated from the fact that populations 
must be discredited. Furthermore, this force macroscopically jeopardizes the biodiversity 
of an ecosystem. However, populations in an ecosystem are discrete integers. Approxi- 
mating these discrete populations by continuous variables inevitably introduces intrinsic 
fluctuations, which turn the evolutionary dynamics stochastic in nature. When external 
noises are introduced, the stochastic process has been shown to be capable of causing mass 
extinction [13] in analogous to the avalanche in the sand piles. How important is this tiny 
difference between discrete and continuous variables when the population size is large? Do 
the intrinsic fluctuations simply introduce small irregularities or will they ever accumulate 
and cause a drastic impact on the biodiversity of the ecosystem? 

To put the discussions on firm ground, we concentrate on the non-transitive rock-paper- 
scissors game [TH-ITU] . known as a paradigm to illustrate the species diversity. When three 
subpopulations interact in this non-transitive way, we expect that each species can invade 
another when its population is rare but becomes vulnerable to the other species when over 
populated. The non-hierarchical competition [201 - I23] gives rise to the endlessly spinning 
wheel of species chasing species and the biodiversity of the ecosystem reaches a stable dy- 
namical balance. This cyclic evolutionary dynamics has been found in plenty of ecosystems 
such as coral reef invertebrates [24] . lizards in the inner Coast Range of California |2S] and 
three strains of colicinogenic Escherichia co/z [15], [26] in Petri dish. Although the oscilla- 
tory solutions for the replicator equations capture the main features, inclusion of mobility [7] 
or/and finite-population effects[22l [23] in the numerical simulations always jeopardizes the 
stable equilibrium and highlight the importance of stochasticity in the evolutionary dynam- 
ics. 

To measure the effects due to the discreteness of the populations, we introduce a bio- 



diversity indicator which is a direct product of all three subpopulations to specific powers 
and remains constant in the continuous replicative evolution. By extensive numerical simula- 
tions, we record how the biodiversity indicator receives random corrections from the intrinsic 
fluctuations. In addition to the irregular deviations at the short-time scale, it is truly re- 
markable that dissipative dynamics emerges as the stochastic processes accumulate and the 
biodiversity indicator thus decays exponentially. Slowly but surely, one species will flrst 
become extinct after a half-time proportional to the population size|27t 128] . Our flndings 
can be elegantly summarized in three steps: discreteness induces fluctuations, fluctuations 
spawn dissipations and dissipative dynamics leads to extinction. 

The subtle connection between fluctuations and dissipations is best exemplifled by a 
damped simple harmonic oscillator moving in a viscous liquidp9]. The microscopic random 
bombardments from the thermal molecules coarse-grain into a macroscopic friction, causing 
the otherwise-conserved mechanical energy to dissipate exponentially. Although our analysis 
is based on the cyclic-competing ecosystem, the emergent dissipative dynamics from the 
intrinsic fluctuations can be readily applied to general biological and ecological systems. It 
can also be generalized to many other practices, such as the Turing [30] model for biological 
patterns [3 Ij or quantitative economic modelling on capital stock [32] and reverse logistics [33] 
where the basic ingredient is respectively the discrete pigment and monetary unit. 

Here we start the exploration on the rock-paper-scissors game and investigate how the dis- 
sipative dynamics arises from the discreteness of the populations. Consider three cyclically 
competing species A, B, C with the stochastic interactions among them, 

A + B -H A + A, 
B + C ^ B + B, 

C + A ^ C + C, (1) 

where kab, kbc, kca are the relative probabilities for the cyclical replacements to occur. For 
convenience, we choose the normalization kab+kbc+kca = 1- Not that the stochastic processes 
preserve the total population A^ = Na + Nb + Nq, where Na, Nb, Nq are the subpopulations 
for each species. In each simulation time step, a pair of individuals is randomly chosen and 
evolves stochastically according to the cyclical competing processes. 

To establish the connection between the stochastic and the replicative approaches, it is 
insightful to derive the effective replicator equations for the population ratios deflned as 
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FIG. 1: Evolution for the cyclically competing species from replicator equations. (A) The cyclic 
competition generates the delicate dynamical balances among the species and leads to the stable 
oscillatory populations for the competing species. (B) The oscillatory population densities trace 
out a closed simplex contour (purple solid line) inside the triangle. For finite population size 
A^ = 9, the changes of the populations ratios are described by stochastic directional hopping on 
the underlying triangular lattice inside the simplex. 

a{t) = NA{t)/N and similarly for b{t) and c{t). As long as the total population A^ is large, 
we can approximate the time derivatives of the population ratios as [(a(t + r)) — {a{t))]/T. 
The detail derivations for the effective replicator equations scratched from the microscopic 
stochastic processes can be found in Supplementary Online Materials (SOM), 

d = a{kabh-Kac) +ia{t), 
b = b{kbcC - kabO) + ^b(t), 
c = c{kcaa - hch) + ic{t), (2) 

where ^a(t) , $,b(t) , ^c(t) are intrinsic noises arisen from the discreteness of populations and 
only become insignificant in the limit of infinite population. Even though these intrinsic 
noises are not correlated at different times, they are not completely independent. Because 
Na + Nb + Nq is strictly conserved at each microscopic evolution step, a global constraint 
holds among these noises, iaif) + ibif) + icif) = 0. 

To visualize the origin of the intrinsic fluctuations, it is convenient to plot the evolution 
trajectory in an appropriate phase space. Let us ignore the noise terms in the replicator 
equations momentarily and study the symmetric case with kab = kbc = kca = 1/3. It is 
straightforward to solve the replicator equations and obtain the oscillatory solutions for the 
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FIG. 2: Comparison between the stochastic evolution of Na — Nf, (color solid lines) with finite 
population size N = 500 and that from the replicative dynamics (grey lines). (A) Starting inside 
the simplex ^3, the stable oscillatory solution from the replicator equations is gradually damped 
out by the intrinsic fluctuations. (B) On the absorbing boundary, the cyclic competition no longer 
reigns and the competing mechanism is the usual prey-predator type. The population evolution in 
this hierarchical case can be solved analytically for different kab = 2,3,4.5,8 (from left to right in 
grey colors). The stochastic evolution for the hierarchical competition gives identical results (green, 
yellow, orange and red curves) to the replicative dynamics. As shown in the inset, all evolution 
trajectories show nice scaling behavior and can be collapsed onto a universal curve when plotted 
against the rescaled time kabt. 



population ratios as shown in Fig. lA. Due to the constraint a + b + c = 1, the allowed phase 
space for the evolution trajectories comprises a standard triangle [TT]. known as the simplex 
5*3 with three absorbing corners where only one species survives. In addition, along the three 
boundaries where one competing species becomes extinct, the evolution dynamics is just the 
usual prey-predator type and not cyclic in nature anymore. The oscillatory solution for the 
replicator equations traces out a closed contour in Fig. IB clockwise inside the simplex 5*3. 
When the size of population is finite, the changes of population ratios occur in steps of stride 
length 1/A^ and form a microscopic triangular lattice inside the simplex. For visual clarity, 
the underlying triangular lattice for a small population A^ = 9 is plotted in Fig. IB. The 
stochastic processes lead to zigzag hopping on the triangular lattice, while the replicator 
dynamics predicts a smooth contour inside the simplex. The deviations between the zigzag 
hopping and the smooth trajectory is the origin of the intrinsic noises. It shall be clear 
that the strength of these noises is of order 1 /N and will only vanish in the limit of infinite 
population. 

One may naively expect that the noises make the evolution wandering around the smooth 
trajectory predicted by the replicative dynamics in a random fashion. If so, as long as the 
trajectory is far from the absorbing boundaries of the simplex, the random wandering can be 
averaged out. We choose an initial set of the population ratios inside the simplex and solve for 
the evolution by both stochastic and replicative approaches. The differences between these 
two approaches are compared in Fig. 2A. The replicator equations deliver the oscillatory 
solutions for the population ratios as expected. In contrast, the stochastic processes seem 
to damp out the indefinite oscillations fairly quickly and lead to extinction with only one 
surviving species. The drastic difference shows that the shot noises provides more than 
mere random drags to the smooth contour, as is predicted by the replicative dynamics and 
appropriate when the initial set happens to fall on the absorbing boundaries. 

To capture the non-trivial effects of the intrinsic noises, we introduce a biodiversity indi- 
cator 

X{t) = a{tf^%tf-c{tf'^\ (3) 

which is constant within the replicative dynamics and serves as a good quantitative mea- 
sure for how the biodiversity enervates during stochastic evolution. Following the standard 
derivations from the fluctuation-dissipation theoremj29]. we arrive at the central result of 




FIG. 3: The biodiversity indicator xi^) of the cychc- competing ecosystem. (a) Several 
single-round evolutions in the numerical simulations are shown in different colors, where 
{kab,khc,kca)=iO -30, 0.35, 0.35) and N = 900 with equal initial subpopulations. The biodiversity 
indicator is constant within the replicative dynamics (purple solid line) but driven to extinction 
X = by random events in the stochastic evolution, (b) After ensemble averaging through 1000 
samples, the biodiversity indicator is found to converge to a smooth exponential decaying function, 
as corroborated by the straight line on the semi- logarithmic plot in the inset. 



the paper. 



d{x) _ {Xit + t)) ~ ixit)) 
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Dia,b,c)(3{x), 



(4) 



dt T N 

where the detail steps can be found in SOM. The factor D{a,b,c) conies from the noise 
correlator for x and is thus positive at all times. The 1/A^ factor is pulled out to emphasize 



that the strength of the noise correlator is inversely proportional to the total population. 
The key player for driving the dissipative dynamics is the asymmetry in the phase space, 

where L{x) denotes the contour length for the particular x- Note that x reaches its maximum 
at the fixed point (a, b, c) = {k^c-, kca, kab) and decreases monotonically when approaching the 
absorbing boundaries. It is clear that L = at the fixed point and increases monotonically to 
L = 3 when approaching the boundaries. As a result, dL/dx < 0, implying f3{x) is positive- 
definite. Basically, the coarse-grained fluctuations generate a dissipative drive toward the 
absorbing boundaries where the phase space is largest. It is rather remarkable that the 
direction of dissipative evolution is dictated by the geometric structure of the simplex S3. 
When one species is removed from the competition, the phase space is reduced to one of the 
absorbing boundaries. The same reasoning can be applied to the hierarchical competition 
between the remaining two species and predicts no preferential drift direction since the one- 
dimensional phase space is symmetrical. The absence of dissipative drive on the absorbing 
boundaries agrees with the numerical simulations to be explained later. 

Numerical simulations are performed to back up our predictions. The ensemble-averaged 
x{t) in Fig. 3 indeed reveals the monotonic decreasing trend. Except for the initial transient 
period, the dissipative dynamics is well captured by the exponential decay. This universal 
exponential decay is due to the simplification of the dissipative equation around small x{t)- 
When close to the absorbing boundaries, the product D{x)P{x) is small and mainly depends 
on X- We can then expand the equation near x = to the linear order, 

d{x) f 1 



dt V Ntsx J 

where I/t^x = d{Dl3)/dx\x=o- As a result, the extinction process is exponential in the long- 
time limit and the extinction time T^x = Nt^x of the ecosystem scales linearly with the 
system size A^. We perform extensive numerical simulations in the parameter space and find 
the exponential form robust. It is remarkable that the time evolution of the biodiversity 
indicator shows interesting scaling behaviors with respect to the population size as indicated 
in Fig. 4. For a chosen set of competing parameters {kab-, kbc, kca), the biodiversity indicators 
for different population sizes A^ can be collapsed onto a universal curve when plotted versus 
the rescaled time t/N. Note that the rescaled extinction time Tex = T^x/N is closely related 
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FIG. 4: The time evolution for the biodiversity indicator x with different competing parameters and 
population sizes. Here we choose three sets of {kab, khc, kca) = (0.30, 0.35, 0.35), (0.15, 0.20, 0.65) and 
(0.10,0.45,0.45), corresponding to the middle (green), corner (red) and boundary (blue) contours. 
The population sizes are N = 90, 180, 360, 540, 720, 900 with equal initial subpopulations. It is 
rather impressive that the biodiversity indicator always decays exponentially. Furthermore, the 
scaling behavior for different population sizes is manifest when plotted versus the rescaled time 
t/N. All curves with the same competing parameters fall onto the universal curves, captured well 
by the exponential decays (dashed lines). 

the location of the evolutionary trajectory. As shown in Fig. 4, for contours deep inside the 
simplex, the dissipation takes a longer time to drain up the biodiversity. On the other hand, 
the dissipation is understandably much quicker when contours are close to the boundaries 
or corners. However, despite of the differences in details, the almost-perfect collapse shows 
the universal exponential decay with 1/A^ scaling. 



The macroscopic dissipation emerges from the microscopic fluctuations is best under- 
stood in the context of fluctuation-dissipation theorem[2H] for general dynamical systems. 
Table I provides a clear similarity of the shared properties between the cyclic competition 
and the damped simple harmonic motion (SHM). However, there are still marked differ- 
ences which make the dissipative dynamics in the flnite-population ecosystem unique. In 
ordinary dynamical systems, the liquid reservoir provides the source of fluctuations and also 



shared properties 


cyclic competition 


damped SHM 


conserved quantity 


biodiversity indicator x 


mechanical energy E 


nature of fluctuations 


intrinsic 
discreteness of populations 


extrinsic 
random colhsions from reservoir 


resulting dissipation 


dt ^^ 


f <o 


eventual outcome 


definite extinction 


come to a stop 



TABLE I: The close analogy between the cyclic-competing ecosystem and the damped simple 
harmonic motion (SHM) in a viscous liquid. 

establishes the asymmetry of the phase space since the degrees of freedom in the reservoir is 
much larger than those of the studied dynamical system. However, in the cyclic-competing 
ecosystem, there is no external reservoir and the fluctuations are intrinsic and originated 
from the discreteness of species populations. In addition, the asymmetry of the phase space 
for the ecosystem is not always guaranteed because there is no external reservoir. By now, 
we have established a new source of dissipations to the non-transitive competitions inside 
the simplex S^. Although these dissipations are absent for the hierarchical competitions on 
the absorbing boundaries, they may be resurrected when the structure of the phase space 
is changed upon the inclusion of other competing mechanisms. An example of this is the 
Lotka-Volterra model in which both selection and reproduction of preys and predators are 
taken into account. The evolution trajectory becomes two dimensional and the phase space 
shares the same asymmetrical structure as for the cyclic-competing ecosystems. We thus 
expect similar dissipative dynamics will occur and bring down the biodiversity. It is rather 
profound that the asymmetrical structure of the phase space pins down the direction for the 
dissipative evolution toward extinction. 

In conclusion, we investigate the intrinsic fluctuations from the discreteness of the popu- 
lations and show how the dissipative dynamics for the biodiversity indicator emerges. Our 
findings pave a different route to address the intrinsic noises beyond the replicator dynamics 
and deepen our understanding for the stability of biodiversity and the evolutionary dynamics 
toward extinction in ecosystems with finite populations. 
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INTRINSIC FLUCTUATIONS 



We would like to derive the intrinsic noises from the discreteness of the populations. For 
instance, the population for species A after each stochastic evolution step changes AA^^ = 
1,0,-1 with probabilities Pi,Pq,Pi respectively. We can separate the changes into the 
average part and the fluctuations. 



ANair) = Rab - Rca + Ur) 



(6) 



where Rab = kabdh and R^a = kcaCa. The simulation time r = 1,2,3,... denotes the gen- 
eration of the stochastic evolution. By construction, the probability distribution for the 
intrinsic noise is 



Rr 



Rr 



Pi = Rab 

Pi = Rca ~> 'Ca = — 1 — {Rab ^^caj, 

Pq = I — Rab — Rca — ^ ^a = — (Rab — Rca)- 



(7) 
(8) 
(9) 



It is easy to check that the average of the noise is indeed zero. The correlation function for 
the noise is straightforward to compute. 



{Ur)Ur')) = Srr' E PiX)ea{X). 

X=1,1,0 

After some algebra, the noise correlator for species A is 



(10) 



{Ur)Ur')) = 5r 



Rab + Rca — {Rab — Re 



Similarly, one can compute the other noise correlation functions for species B and C, 



UrMr')) = 5r- 
{Ur)Ur')) = Sr 



Rbc + Rab ~ {Rbc — R, 



ab 



Rca + Rbc — {Rca — R 



be 



(11) 

(12) 
(13) 



Suppose the population size is large enough that the finite difference can be approximated 
by the continuous differential, 

da Aa ANa 



dt At At 



Rab — Rca + (,a{t), 



(14) 
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where Aa = ANa/N and At = At/N = 1/N. Therefore, the rephcator equation is recovered 
with an intrinsic noise ^a due to the discreteness of the population. Note that the strength 
of the intrinsic noises exhibits a 1/A^ factor in the continuous hniit, 

1 



{Ut)Ut')) ^ 



Rab + Rca — {Rab — R 



caj 



5{t-t'). (15) 



The rephcator equations with intrinsic noises for species B and C can be derived hkewise. 

DISSIPATIVE DYNAMICS 

In the following, we shall borrow the lesson of fluctuation-dissipation theorem in 
physics [29] to show the emergence of the dissipative dynamics for the biodiversity indicator. 
Since the biodiversity indicator x(t) is a constant of motion for the replicator equations, 
the dynamics is driven purely by the corresponding noise ^x(^)- Integrating the dynamical 
equation for x{i)i the finite difference after ensemble average is 

{X{t + r))-{x{t)) = [^^dt'{Ut'))- (16) 

It is important to emphasize that r is much smaller than the macroscopic time scale we are 
interested in, but much larger than the microscopic time scale for evolution steps, 

T* <t:T<t: At. (17) 

Since t ^ t*, the ensemble average makes sense at time t' with a modified probability 
distribution P{t') which differs from the original one P{t) = Pq. 

Assuming the probability distribution function P{t') is proportional to the corresponding 
phase space L{x + Ax), 

Pit') L{x + Ax) ^^ _^A^ 



P(t) L{X) 



e-^^>^. (18) 



t' 



The finite difference in the biodiversity indicator is Ax = /^ ^^{t")dt" and the function 
/3(x) = —d{\nL)/dx > because L{x) decreases as x increases. Putting all pieces together, 
we obtain the important identity 

iUt')) = {Ut')e-'^'')o (19) 
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The ensemble average can be computed by expanding the exponential, 

i^xit')) ^ (ex(0(l-/3Ax))o = -/3(x)(ex(OAx)o 

= -/3(x) I' dt"{i^{t')Ut")), = -l/3(x)D(a, 6, c), (20) 

where both P{x) and D{a,b,c) are positive definite. Finally, the time evolution of the 
biodiversity indicator is 

d{x) _ {X{t + r)) - ixit)) 1 



dt T N 



/3{x)D{a,b,c)<0. (21) 



It is quite remarkable that we prove that the ecosystem tends to evolute toward extinction 
after coarse-graining the intrinsic fluctuations from the discreteness of the populations. 
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